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The limits of applicability of the Lang-Kobayashi (LK) model for a semiconductor laser with op- 
tical feedback are analyzed. The model equations, equipped with realistic values of the parameters, 
are investigated below solitary laser threshold where Low Frequency Fluctuations (LFF) are usually 
observed. The numerical findings are compared with experimental data obtained for the selected 
polarization mode from a Vertical Cavity Surface Laser (VCSEL) subject to polarization selective 
external feedback. The comparison reveals the bounds within which the dynamics of the LK can be 
considered as realistic. In particular, it clearly demonstrates that the deterministic LK, for realistic 
values of the linewidth enhancement factor a, reproduces the LFF only as a transient dynamics 
towards one of the stationary modes with maximal gain. A reasonable reproduction of real data 
from VCSEL can be obtained only by considering noisy LK or alternatively deterministic LK for 
extremely high a-values. 

PACS numbers: 42.65.Sf,42.55.Px,05.45.-a,42.60.Mi 



I. INTRODUCTION 



The dynamics of semiconductor lasers with optical feedback is studied both experimentally and theoretically since 
almost 30 years (0], for a review see e.g. Q ). The interest for such configuration, commonly encountered in many 
applications (e.g, communication in optical fibers, optical data storage, sensing etc) arises from the rich phenomenology 
observed, ranging from multistability, bursting, intermittency, irregular and rare drops of the intensity (low frequency 
fluctuations (LFF)) and transition to developed chaos (coherence collapse (CC)). A complete understanding of the 
physical mechanisms at the basis of such complex behavior is, however, still lacking. In particular, the origin of the 
LFF regime is under debate since the very first observations and yet this puzzling problem has not been solved. Their 
origin was ascribed to stochastic effects 0, or to deterministic but chaotic dynamics |f| , and more recently even to 
the interplay between regular periodic and quasi-periodic solutions [||. The LFF dynamics has been investigated by 
using several type of emitters, mainly edge-emitting: ranging from longitudinal multimode 0, IE U EH El C3 E3 to 
single-mode DFB |l4| semiconductor lasers. 

From the experimental point of view, a complete characterization of the LFF dynamics is quite difficult, because 
of the very different time-scales involved |l5|. Indeed, fast oscillations on the 10-ps range have been observed with 
streak camera measurements [l6l[l7j . representing the fundamental scale on which the system evolves. On the other 
hand, the duration of such fast pulsing regime between LFF events can be as long as hundreds of nanoseconds or even 
microseconds. In literature, the LFF dynamics has been experimentally characterized in several manners, starting 
from a relatively simple statistical analysis of the time separation T between LFF0, 0, Il9ll (relating the average 
(T) between LFF with the pump current) to Hurst exponents for the laser phase dynamics [201 ]. 

A widely used theoretical description of the system is the Lang-Kobayashi model [21] , introduced in 1980 in the effort 
to provide a simplified but effective analysis of an edge-emitting semiconductor laser optically coupled with a distant 
reflector. In the model, both the multiple reflections from the mirror (low coupling) and the possible multimodal 
structure of the laser were neglected. The opportunity to include such effects has been discussed in several papers 
0, B S E3i E3 7 but the model still remains presented as the standard theoretical approach to the system. While 
most of the phenomenology observed in the different experiments is grab by the model, quite often a more precise 
or quantitative comparison is obtained at the expense of a choice of parameters far from those actually measured or 
even not physically plausible. 

Recently, a new configuration has been proposed and studied, based on a Vertical Cavity Surface Laser (VCSEL) 
with a polarized optical feedback ^5|. Such a laser is longitudinal single-mode (due to the very short cavity) but 
may support different, high order transverse modes for strong enough pumping current (see e.g.[24|). The simmetry 
of the cavity allows also for the possible laser action on two different, linear polarizations selected by the crystal 
axis. The dynamics of the VCSEL with isotropical optical feedback has been examined experimentally in H^H^] and 
theoretically in [2(|, while the role played by polarized optical feedback has been discussed in [2?ll2^.l29| . In particular, 
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the setup used in jl5| employed a polarizer in the feedback arm, in order to couple back only the radiation of one 
polarizations; moreover, a suitable range of pump current was chosen, to assure single transverse mode behavior. In 
such configuration, the appearance of LFF was reported and characterized. The possibility to control the role of the 
laser modes in the dynamics in this setup allows for a consistent description via the LK model and therefore for an 
effective test of its predictions, at variance with similar setup where instead the polarization selection were not used 

MM 

Our aim in the present paper is to clarify the origin of the LFF dynamics by comparing numerical results obtained by 
solving the LK equations with experimental measurements done on a VCSEL. In particular, the parameters employed 
for the integration of the LK rate equations have been carefully derived from the analysis of the same VCSEL employed 
for the experiments [32l| . 

In Sec. II we describe our experimental setup, reporting the main phenomenology observed in the range of variation 
of the more relevant parameters of the system, namely, the pump current and the phase of feedback. The LK model 
is introduced and commented in Sec. Ill, together with the numerical methods employed for its integration and the 
choice of the parameter values derived from the experiment. In Sec. IV the properties of the stationary solutions 
are discussed, while in Sec. V a careful characterization of the deterministic model is given, detailing the transient 
phenomena and the Lyapunov analysis. In Sec. VI the effect of noise is introduced and analyzed, discussing also 
its possible importance in the experiment. A detailed comparison of the numerical results with the experimental 
measurements is given in Sec. VII, with particular regard to the distribution of the intensity and of the inter-event 
times for different parameter choice, including the alpha-factor and the acquisition bandwidths. Finally, we draw our 
conclusions in Sec. VIII. 



II. EXPERIMENTAL SETUP AND SETTINGS 




I (mA) Time (ns) 

FIG. 1: (a) Average power output versus the input current for the solitary laser (dots) and for the laser with feedback (stars). 
(Color online) (b) Polarization modes: power outputs as a function of time for the VCSEL with feedback at I = 2.64 mA. 
Upper trace (black): main polarization; lower trace (red): secondary polarization. 

The experimental measurements are performed using a VCSEL semiconductor laser with moderate polarized optical 
feedback. In particular, our analysis is limited to a regime of pumping below the solitary threshold I t h ~ 2.76 mA, 
where the VCSEL emits light on a single linearly polarized transverse mode. Longitudinal modes are not allowed by 
the cavity, and other transverse modes are not present up to currents ~ 6.5 mA. The solitary laser emission remains 
well polarized up to roughly the same current (see Fig. ^(b)). 

The technical details of the source are the following. The laser is an air-post VCSEL made by CSEM |30j , operating 
around 770 nm. The mesa diameter is 9.4 //m, a ring contact defines the ouput window with a diameter of 5 /im, and 
the active medium is composed of three 8 nm quantum wells. The temperature of the laser case is stabilized within 
1 mK, the pump current is controlled by a home made battery operated power supply whose current noise is below 
40 pA/Hz -1 / 2 in the frequency range from 1kHz to 3 MHz. 

The feedback is applied on the polarization direction of the solitary laser emission. The external cavity includes 
collimation optics, two polarizers, a variable attenuator, and the feedback mirrors that is mounted on a piezoelectric 
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transducer at about 50 cm from the laser. The output radiation, after optical isolators, is detected by an avalanche 
photodiode with a bandwidth of about 2 GHz, whose signal, sometimes after low- pass filtering, is recorded by a 4 GHz 
bandwidth digital scope. More details on the experiment can be found in Refs |15ll3ll|. 

Optical feedback results in a reduced threshold JT^ ~ 2.42 mA, as shown in Fig. ^( a )- We have examined the 
dynamic behavior of the output intensity for various pump currents, both above and below I t h, with a particular 
attention to possible effects of the feedback phase A<f> which is varied by acting on the external mirror piezoelectric 
transducer. 

As shown in Fig. [21 we can identify several regimes: (I) at / < Ith one observes a single mode LFF dynamics (i.e. the 
main polarization exhibits LFFs, while the secondary polarization remains off, see Fig. |I](b)); (II) Ith < I < 3.5 mA: 
in this regime the LFF dynamics of the main polarization is accompanied by a synchronized spiking behavior in 
the secondary polarization (coupled modes LFF). This regime was analyzed in Ref. fl5| and, in more details, in 
Ref. |3lJ j where it is proposed that the dynamics of the secondary polarization is driven by the main polarization, 
whose behavior is not influenced by the orthogonal polarization mode. For larger pump current one begins to observe 
Coherence Collapse (regime (III) in Fig. EJl and moreover the feedback phase begins to play a fundamental role. In 
particular, for increasing I in larger and larger portion of the phase interval the laser is stationary (regime (IV) in 
Fig. |2J). This phenomenon is not yet understood and it will be subject of a future analysis [27j|. 

Anyway, in the present paper we will limit our analysis to regime (I) (i.e. for / < Ith)-, where the VCSEL has a 
single-mode LFF dynamics and the phase delay of the feedback does not play any role. We remark that this statement 
is suggested by the experiment, where the phase stability would be enough to discriminate such effects, as shown in 
the analysis of regimes (III- IV). 




I(mA) 

FIG. 2: Phase diagram of the VCSEL with feedback: phase of the feedback A<f) as a function of the pump current /. The roman 
numbers denote regions of different dynamical regimes: (I) single-mode LFF, (II) two-mode LFF, (III) coherence collapse and 
(IV) stable emission. The vertical dashed line indicates (from low to high current) successive current thresholds: the reduced 
one IXht the solitary laser one Ith and the current value separating LFF from coherence collapse. The dots represent the 
maximal A<f) for which the laser emission remains stable, while the solid line is a guide for eyes to distinguish regions (III) from 
(IV). 



III. NUMERICAL MODEL AND METHODS 

The dynamics of the VCSEL for / < I t h is a purely single mode dynamics and therefore we expect that it could 
be reproduced by employing the Lang-Kobayashi (LK) [21j rate equations for the complex field E(t) and the carrier 
density n(t). In order to achieve an accurate and reliable comparison of the numerical results with the experimental 
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ones we will employ for our simulations the laser working parameters reported in TableQ] These parameters have been 
determined via a series of suitable experiments for exactly the same VCSEL employed to obtain the measurements 
examine in this paper • The only lacking parameter that was not possible to obtain in the previous characterization 
is the feedback strength, which is however determined from the threshold reduction. 

In this article we use the rate equations derived in Ref. [32( f° r the single-mode solitary laser, modified to include 
the feedback. Defining the deviation An of the carrier density from transparency normalized to have unitary value 
at threshold, i.e. 

An = -pi , (1) 

n t h - 1 

the following form is obtained [3^ |: 



r n An = -An+l + r)(jjb- 1) - An\E 
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E = —^{An-l]E+—c- WT E(t~T) + ^R £,{t) (2) 

where £(t) = £ij(t) + i£i(t) is a complex Gaussian noise term with zero mean and correlation given by < 
£fl(i)£-R(0) >= < £r(*)£r(0) > = <K*) an d < £r(£)£t(0) >= 0. The noise variance R = {n/n t h) 2 R S p represents a 
multiplicative noise term proportional to the square of the reduced carrier density n/n t h {nth being the threshold 
carrier density) and to the variance of the spontaneous emission noise R sp . The parameter fi = I j I t h is the pump 
current rescaled to unity at threshold, t„ and t p are carrier and photon lifetimes, respectively, r is the delay (or 
external roundtrip time), a is the linewidth enhancement factor and r/ the reduced gain (for the exact definitions of 
these quantities in terms of the laser parameters and for the approximations employed to derive J2J see Ref. 

By reexpressing the time scale in terms of the photon life-time (r p ) the equations assume the usual form for the 
LK model and read as 

TAn = -An+p~ An\E\ 2 

E = ^^[An-l]E + kc-^ T E{t-T)+VRm (3) 

where T = T n jr v , p = 1 + r;(u - 1) and R = (n/n th ) 2 R sp x t p . 

The complex field can be expressed as E(t) = A(t) + iB(t) = p(t) expiip(t) and the equations can be rewritten as: 



A(t) = 

B(t) = 



^ (B(t) - aA(t)) + k[B(t - t) cos0 - A(t - r) sin</>] + y/R/2^{t) 



An(t) - 1 
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TAn{t) = p- An(t) [l + A 2 (t) + B 2 (t)] (4) 
where cj> = lot. 

Moreover, by assuming that n ~ n t h the noise terms become additive with an adimensional variance R — 2.76 x 10 -3 , 
the other quantities entering in expressed in r p units, are t = 302.5, uj x t = 8.743 x 10 6 and T = 30.8333, the 
numerical values have been obtained by employing the parameters values in Table [fl 

In order to reproduce the power-current response curve for two different experimental data sets we have chosen 
feedback strengths k = 0.25 and 0.35, while typically we considered pump currents and linewidth enhancement factors 
in the range 0.9 < fi < 1.20, and 3 < a < 5, respectively. 

The deterministic equations Q with R = have been integrated by employing the method introduced by Farmer 
in 1982 |33| equipped with a standard 4th order Runge-Kutta scheme, while for integrating the equations with 
the stochastic terms we have employed an Heun integration scheme [34]] . The simulations have been performed by 
integrating the field variables with time steps of duration At = r/ (N — 1), with N = 1, 000 — 10, 000. 

The dynamical properties of the system can be estimated in terms of the associated Lyapunov spectrum, that fully 
characterizes the linear instabilities of infinitesimal perturbations of the reference system. By following the approach 
reported in |33| . we have estimated the Lyapunov spectrum {A^} (k = 1, . . . , 2N + 1) by integrating the linearized 
dynamics associated to eqs Q in the tangent space and by performing periodic Gram-Schmidt ortho-normalizations 
according to the method reported in |35| . The Lyapunov eigenvalues are real numbers ordered from the the largest 
to the smallest, a positive maximal Lyapunov Ai is a indication that the dynamics of the system is chaotic. Moreover, 
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from the knowledge of the Lyapunov spectrum it is possible to obtain an estimation of the number of degrees of 
freedom actively involved in the chaotic dynamics in terms of the Kaplan- Yorke dimension [3(j : 



D 



KY 



|A i+ i| 



where j is the maximal index for which 53i=i > 0. 



(5) 



TABLE I: Experimental values of the parameters entering in the model © (from j^j) 



Description 



Symbol 



Value 
3.2 ±0.1 
12 ± 1 ps 
0.37 ±0.02 ns 
3.63 ns 

(2.3 ±0.6) x 10" 4 ps' 1 
5.8 ±0.6 



Linewidth enhancement factor 

Photon lifetime in the cavity 

Carrier lifetime 

External roundtrip time 

Variance of the spontaneous emission noise 

Reduced gain 



Tn 

T 



lis 



IV. STATIONARY SOLUTIONS 



A first characterization of the phase space of the LK system can be achieved by individuating the corresponding 
stationary solutions and by analyzing their stability properties. The stationary solutions of the above set of equations 
can be found by setting p — Ah = and "0 = f2, i.e. by looking for solutions of the form 

E s (t) = p S e lQt and An(t) = An s . (6) 

The stationary solutions can be parametrized in terms of the variable 9 = <f> + fir, and by setting J = (p — l)/2 they 
can be expressed as follows 37] 

X S = = -fccosW (7) 



and 

n = — - = -k\/l ± a 2 sin(0 ± ) (9) 

T 

where 9q = atan(a) and — n < 9 < ir. 

The stationary solutions are usually termed external cavity modes (ECMs) and correspond to stationary lasing 
states, in the limit of long delay they are densely covering the following ellipse: 

k 2 = X 2 S ± (n - aX s ) 2 . (10) 

The stability properties of these solutions can be obtained by estimating the associated Floquet spectrum {A„}, 
these eigenvalues are typically complex and their number is infinite, due to the delay term present in the LK equations. 
However, the stability properties of the ECMs are determined by the eigenvalues with the largest real part, in particular 
by the maximal one A M = (Re A M ,Im A M ). These can be easily determined by solving the characteristic equation 
obtained by linearizing around a certain ECM: 

Z 2 [A n + e(l + p 2 s )] ± 2Z„[A„ cos(<9)(A„ ± e(l ± p%) + ep%{l - 2k cos(0))(cos(0) - a sin(0))/2] 
+ A 2 n [An ± e(l ± p 2 s )] ± A n ep|(l - 2fccos(0)) = a n Z 2 n ± 2b n Z n + c n = . (11) 
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fc(l 



-A„T 



) and e = l/T. The pseudo-continous spectrum [3^] can be obtained by solving eq. 



where Z n 

in terms of Z n and by considering A„ as a parameter. In this case the solutions of the corresponding second order 
equation are 



A. 



■log 



i2im 



(12) 



and by self-consistently solve the above equation one can find all the eigenvalues, that are arrangend in two branches. 
However, the spectrum contains also isolated eigenvalues, that can be obtained by solving directly Eq. in terms 
of A n and by considering Z n as a parameter. In this case the equation is a cubic and by solving self-consistely its 
expression one finds (up to) three distinct eigenvalues. While the pseudo-continuos spectrum emerges due to the 
presence of the delay in the system, the isolated eigenvalues originate from those characterizing the three dimensional 
single laser rate equations in absence of the delay, i.e., the Eq. © with r = |3S[ . 

Typically, depending on their linear stability properties ECM are divided into modes and antimodes [39l | . Antimodes 
are characterized by a positive real eigenvalue and are therefore unstable. For the modes instead the maximal real 
eigenvalue is zero and they are unstable whenever the associated spectrum crosses the imaginary axis. Various types 
of instabilities can be observed for these delayed systems and they can be classified for analogy with the spatially 
extended systems as for example modulational-type or Turing-type instabilities |3?l l38|. Examples of the unstable 
branch for the spectra {A^} associated to Turing-type and modulationali-type instabilities are reported in Fig. |3 for 
a classification of the possible instabilities of equilibria for delay-differential equations see (3^ . 




0.04 



Re A (ns 



FIG. 3: (Color online) Unstable branch A associated to specific modes of the LK. The mode 9 — —0.32 for a — 3.20 reveals 
a sort of modulational instability (red squares), while the mode 6 = —0.08 for a = 3.22 shows a Turing-like instability (blue 
circles). The results refer to /i = 0.93. Due to the symmetry of the spectra only the part corresponding to ImA > is displayed. 



In particular, the modes correspond to 0- values located within the interval [62 '■ #1], where 0\ ~ atan(l/cv) and 
62 ~ — 7T + atan(l/a). Moreover, stationary solutions are acceptable only if they correspond to positive intensities 
Ps ^ 0; i- e - ^ they are situated within the interval [— Or : Or] with Or = acos(— J/k). In the range of parameters 
examined in the present paper the number of stable modes (SMs) is always between two and four and they are located 
in a narrow interval of located around zero (i.e. around the so called Maximum Gain Mode (MGM)). 

As we will report in the following section, by integrating the deterministic version of the model Q for moderate 
a- and //-values, we always observe a relaxation towards one of the SMs. In particular for 3 < a < 5 the dynamics 
seem to relax always towards one of the SMs located in proximity of the MGM (corresponding to = 0) (see Fig. @}. 
It should be noticed that the MGM is a solution of the system only for appropriate choices of <j>. Two examples of 
typical ECMs are reported in Fig. 0|for the present system for a = 3.2 and 5 and fi — 0.93, in both cases two stable 
attracting modes have been identified. Recently, a similar coexistence of two stable solutions located in proximity of 
the MGM has been reported experimentally for an edge-emitter laser with a low level of optical feedback [40( 

As reported in [^3 the stability properties of the MGM do not depend on a, therefore it is reasonable to expect 
that some SM will be always present in a narrow window around for any chosen linewidth enhancement factor [39| 
and that they will coexist with the chaotic dynamics, as observed experimentally in |4l|. 
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FIG. 4: Intensities versus for ECMs of eq. IgJ with fi = 0.93 and for a = 3.2 (a) and 5.0 (b). The solid line indicates 9\. 
The ECMs with 9 < 6i are modes, while the others are antimodes. The empty circles indicates the SMs towards which we 
numerically observed relaxation of the system by considering up to 100 different random initial conditions. 



V. DETERMINISTIC DYNAMICS 



An open problem concerning the deterministic LK equations is if and for which range of parameters these equations 
faithfully reproduce the experimentally observed dynamics. Particular interest is usually focused on the a-value to 
be employed to obtain a realistic behaviour of the model. 



A. LFFs as a transient phenomenon 

To answer this question we consider the deterministic verson of the model Q (i.e. by assuming R = 0) equipped 
with the parameter values deduced by the experimental data j3^] apart for /j, and a that will be varied. In particular, 
we would like to understand if the system shows a LFF behaviour and if such dynamics is statistically stationary or 
not. In order to verify it, we initialize randomly the amplitude p(t — 0), the phase (j>(t = 0) and the excess carrier 
density An(t = 0), then we follow the dynamics and we examine the evolution of the field intensity p 2 (t). If the 
dynamics end up on a stationary solution we register the time T s necessary to reach it and we average this time over 
many (M) different initial conditions. In order to measure T S1 we estimate the time needed to the standard deviation 
of the intensity (evaluated over sub- windows of time duration t w ) to decrease below a chosen threshold T. The average 
times < T s > are reported in Fig. 

The main result shown in Fig. [3] is that the system after a transitory phase (shorter or longer) settles down to 
a SM and that the duration of the transient increases for increasing a or /i-values. Preliminary indications in this 
direction have been previously reported in [l4|. These results clearly indicate that the LFF dynamics is just a transient 
phenomenon for the deterministic LK equation for commonly employed a-values (i.e., for a ~ 3.0 — 3.5). For each 
fixed n it seems that < T s > diverge above some "critical" a-value, however we cannot exclude that the system 
will also in this case finally converge to a SM. At least we can consider the reported critical values as lower bounds 
above which LFF could occur as a stationary phenomenon and not as a transitory state. With the chosen numerical 
accuracy and due to the available computational resources, for any pratical purpose a transient longer than 0.1 — 1 s 
can be considered as an infinite time. 

An interesting feature is that the variability of < T s > with a is quite wild and it reflects the stability of the modes 
in proximity of the MGM, as shown in Fig. However the average values < T s > (averaged also over a-intervals of 
width 0.1) still indicates a clear trend of the transient times to increase with a. 

We expect that the time scale associated to the convergence to a stable mode is ruled by the eigenvalue with the 
maximal real part (apart the eigenvalue zero, that is always present due to the phase invariance of the LK equations). 
In particular we expect that the intensity of the signal will converge towards the stable mode as 

p 2 (t) ~ P 2 (0)e 2Rc A " t cos(2Tm A M t) + p| 

where A M is the eigenvalue with maximal (non zero) real part associated to the considered stable solution. For small 
\x we have observed that the dynamics can collapse to different stable solutions (typically, from 1 to 3). By estimating 
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FIG. 5: (Color online) Logarithm of the transient times < T s > as a function of a for various ^i-values: namely, /i = 0.93 
(black filled circles), 0.95 (red empty triangles), 0.97 (blue asteriskes) and 0.99 (green empty squares). The data refers to 
k = 0.25, T = 10 -5 , t w — lOOOr = 3.630/^s and M ~ 10 — 100, and they have been obtained by employing an integration time 
step At — 3.63ps. The time scale of the figure is to = 1 ns. The bars reported for each measured value indicate the range 
of variability of < T s > (within a a-interval of amplitude 0.1) due to its finer structure shown in the inset. In the inset are 
displayed the transient times for fj, = 0.93 reported at a higher resolution in a, namely for a resolution of 0.02. 



the probability P m = N m /M to end up in one of these states (being N m the number of initial conditions converging to 
the m-mode) and by indicating the corresponding eigenvalue with maximal real part as A M (m), a reasonable estimate 
of < T s > is given by 

T est = — E |ReAM( m )| ' (13) 

where T is the employed threshold. As it can be seen in Fig. HO the estimation is quite good and the periodicity of 
the two quantities is identical in the examined range of a-values and for fi = 0.93. The expression always gives 
a good estimate of < T s > for a < 4 and for [i < 1, but the agreement worsens for increasing a-values. 

However, even a more rough estimate is capable to give a reasonable approximation of < T s > and in particular to 
capture its periodicity. This estimate is simply given by 

_ lnT 1 

1_ 2 |ReA M | ' ( ] 

where A M is the eigenvalue with maximal non zero real part associated to the stable mode with higher 9 (i.e. the 
first stable mode located in proximity of the anti-mode boundary). 

From the present analysis it emerges that the stable modes play a relevant role for the (transient) dynamics of the 
deterministic LK equations at least for a < 4. Moreover the observed strongly fluctuating behaviour of < T s > as 
a function of a indicates that the choice of this parameter is quite critical, since a small variation can lead to an 
increase of an order of magnitude of the transient time[44| 



B. Lyapunov Analysis 

We have characterized the transient dynamics preceeding the collapse in the stationary state in terms of the 
maximal Lyapunov Ai and of the associated Kaplan- Yorke (or Lyapunov) dimension Dky- In particular, these 
quantities reported in Fig. [7] have been estimated by integrating the linearized dynamics for a sufficiently long time 
period Ti„ t and by averaging over M different initial realizations. 

It is clear from the figures that the average maximal Lyapunov exponent < Ai > is dcfinctely not zero for all 
the considered situations and that it increases (almost steadily) with the parameter a as well as with the pump 
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FIG. 6: (Color online) Transient times < T s > (black filled circles), T es t (red empty triangles) and T\ (blue empty square) 
in nanoseconds as a function of a. The results have been obtained by examining the decay of p 2 (t) with F — 10 -5 and 
t w = lOOOr = 3.630/xs and M ~ 100 — 500. The employed integration time step is At = 0.363 ps and the data refer to k = 0.25 
and = 0.93. 



parameter. Moreover, also this indicator reflects the stability properties of the SMs located in proximity of the MGM, 
by exhibiting large oscillations as a function of a, as shown in the inset of[7fa). 

The values of < Dky > reported in|7fb) clearly indicate that the system cannot be described as low dimensional, 
even during the transient and even below the solitary laser threshold. As a matter of fact the number of active degrees 
of freedom ranges between 10 and 50. It should be noticed that < Dky > is determined by the instability properties 
not only of anti-modes and but also of modes that have bifurcated (via Hopf instabilities) becoming unstable for 
increasing a. 
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FIG. 7: (Color online) (a) Average maximal Lyapunov exponents < Ai > as a function of a and for various ^-values below 
threshold: namely, \x = 0.93 (black filled circles), 0.95 (red empty triangles), 0.97 (blue asterisks) and 0.99 (green empty 
squares). The bars reported for each measured value indicates the range of variability of < Ai > due to its finer structure as 
measured within a a-interval of width 0.1. In the inset the data for < Ai > are reported for a higher resolution in a (namely 
0.02) for jj, = 0.93. (b) Average Kaplan- Yorke dimensions < Dry > as a function of \i for a — 4 (black filled circles) and 5 
(red empty squares). All the data refer to k = 0.25, for the < Ai > estimation At = 0.363 ps, M — 500, and Tint = 3.63 ms 
while for the < Dky > evaluation At = 3.63 ps, M — 20 and T; nt = 0.14 ms. 



VI. NOISY DYNAMICS 



We have examined the dynamics @ for increasing level of noise, namely for 10 6 < R < 10 2 . Also in this case, 
for a = 3.3 and for noise levels smaller than 10~ 3 the LFF dynamics only occurs during a transient. However for 
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increasing R values we observed a transition to sustained LFF and the transition region was characterized by an 
intermittent behaviour. These behaviours are exemplified in Fig. |H1 for a = 3.3 and fj, — 0.97. As shown in Fig. 
IHfa) the orbit spends long times in proximity of one of the SM and then, due to noise fluctuations, escapes from the 
attraction basin associated to the stable solution and exhibit LFFs before being newly reattracted by the SM. This 
intermittent dynamics can be interpreted as an activated escape process induced by noise fluctuations, and therefore 
the average residence time < T res > in the attraction basin of the SM can be expressed in the following way 

< T res >cx exp [W/R] (15) 

where W represent a barrier that the orbit should overcome in order to escape from the SM valley. As shown in Fig. 
[5] the process can be indeed interpreted in terms of the Kramers expression (|15|l for 2 x 10~ 4 < R < 7 x 10~ 4 . It 
means that for R < W one should expect an intermittent behaviour, while for R > W the dynamics of the orbit will 
be essentially diffusive, since the noise fluctuations are sufficient to drive the orbit always out of the SM valley. These 
indications suggest that in order to observe a "non transient" or "non intermittent" LFF dynamics the amount of 
noise present in the system should be larger than W. Also in [14( it has been clearly stated that the LK-equations 
with parameters tuned to reproduce the dynamics of a DFB laser with a — 3.4 can give rise to stationary LFF only 
in presence of noise. 

It is important to remark that for a — 3.3 the experimentally measured variance of the noise R = 2.76 x 10~ 3 is 
above the barrier W = 1.87 x 10~ 3 found from the fit of the numerical < T res > with expression Ijl5|l . performed in 
the small noise range (see Fig. [Sj. Moreover in the experiments we never observed relaxation of the dynamics towards 
a SM. 
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FIG. 8: Intensity of the field p 2 as a function of the time for the noisy LK equations. The data have been filtered with a low-pass 
filter at 80 MHz and refer to /i = 0.97, a = 3.3. The results reported in (a) correspond to a noise variance R = 3 x 10" 4 , while 
those in (b) are relative to R = 3 x 10~ 3 . In (a) a typical time of residence T res around one of the SMs is indicated. The inset 
in (b) is an enlargement of the actual dynamics. 



A. Lyapunov Analysis 

Also in the noisy case we have examined the degree of chaoticity in the system by estimating the maximal Lyapunov 
exponent along noisy orbits of the system. The role of noise is fundamental in destabilizing the dynamics of the system 
and in rendering the asymptotic dynamics chaotic. 

In particular, as shown in Fig. HOf a - ) we observe that at /i = 0.93 and k = 0.25 the deterministic dynamics (R = 0) 
is asymptotically stable in the range a < 4.4, while the noisy dynamics becomes more and more chaotic for increasing 
R. For the value R = 3.3 x 10~ 3 , close to the experimental one, the dynamics is completely destabilized in the whole 
examined range 2.5 < a < 4.5, while for smaller i?-values the range of destabilization is reduced. These results confirm 
the role of the noise in rendering the LFF an asymptotic phenomenon. Moreover the maximal Lyapunov increeases 
steadily with a at R = 3.3 x 10~ 3 . Similar findings apply in the case /i = 0.97 (see Fig. HOf blh 

The wild oscillations in the < Ai > values observable at level of noise R < 10~ 4 reflect the stability properties of 
the SMs attracting the asymptotic dynamics. 
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FIG. 9: Average residence times in the SM as a function of the inverse of the variance of the additive noise to the LK equations. 
The dashed line is a exponential fit oc exp [W/-R] to the numerical data, the fitted exponential slope is W = 0.00187. The data 
refer to [i = 0.97 and a — 3.3. 
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FIG. 10: (color online) Maximal Lyapunov exponents Ai as a function of a for /i = 0.93 (a) and /i = 0.97 (b) and for various 
noise amplitude R: namely, the data for R = 3 x 10 -5 are indicated by green asterisxs, those for R = 3 x 10~ 4 by blue filled 
triangles, and the ones corresponding to R = 3 x 10 -3 by red filled circles. The values estimated during the transient dynamics 
in absence of noise are indicated by black empty squares, while the asymptotic values for R — by black filled squares. All 
the data refer to k — 0.25, for the estimation of Ai in the noisy case one orbit has been followed for a time t — 3.63 ms with 
time step At = 0.363 ps, while in the deterministic case the asymptotic results have also been averaged over M = 10 different 
initial conditions. For details on the estimation of the transient Lyapunov exponents see the previous section IV Bl 



VII. COMPARISON BETWEEN EXPERIMENTAL AND NUMERICAL DATA 



This Section will be devoted to a detailed comparison of numerical versus experimental results with the aim to 
clarify if the deterministic or noisy LK eqs are indeed able to reproduce the experimental findings. 



12 



P(P ) 





2 V : 









0.5 



1 



1.5 



2 



2 



P 



FIG. 11: Field intensities distributions P(p 2 ) for the experimental signal filtered at 200 MHz. The experimental intensities 
have been arbitrarly rescaled to match the corresponding average intensities obtained from the simulation of the noisy LK eqs 
at a = 3.3 and k = 0.35 with noise variance R = 3.3 x 10~ 3 . The data refer from left to rigth to I = 2.48 mA, 2.50, 2.54, 2.58, 
2.64, 2.70, and 2.75. Since Rh ~ 2.76 these data correspond to 0.9 < p < 1.0. The first distribution has a large contribution 
from the Gaussian electronic noise, which also explains the negative p 2 values. 



As a first indicator we have considered the distribution of the field intensities P(p 2 ), in particular in order to match 
the experimental findings we consider the signal filtered at 200 MHz. 

The experimental results for the probability distribution functions (PDFs) of the field intensities p 2 are reported in 
Fig. ^] for various currents below the solitary threshold value. It should be noticed that the amplitudes p 2 have been 
rescaled in order to match the corresponding numerical values for the noisy LK eqs with a — 3.3 and noise variance 
R = 3.3 x 10~ 3 , but that no arbitrarly shift have been applied to the data. 

A peculiar characteristic of these data is that for increasing pump current the PDFs become more and more 
asymmetric, revealing a peak at large intensities that shifts towards higher and higher p 2 values and a sort of plateau 
at smaller intensities. 

The corresponding PDFs are reported in Figs. 1121 andl*Hflfor data obtained from the integration of the noisy and 
deterministic LK eqs, respectively. A better agreement between numerical and experimental findings is found for the 
noisy dynamics with a ~ 3.3 — 4.0 and k = 0.35, with a noise variance similar to the experimental one (namely, 
R = 3 x 10~ 3 ). For the deterministic case (reported in Fig. U^for a = 5.0 and k = 0.35) a non zero tail at p ~ is 
observed even for p ~ TO, contrary to what observed for the experimental data. 

These results indicate that it is necessary to include the noise in the LK eqs to obtain a reasonable agreement with 
the experiment, at least at the level of the intensities PDFs. However, the numerical data seem unable to reproduce 
the narrow peak present in the experimental ones at large intensities and for p — ► 1. In the next sub-section a further 
comparison will be performed to validate these preliminary indications. 



We will first compare the experimental and numerical measurements of the average times between two consecutive 
drops of the field intensities < Tlff >• They have been evaluated in two (consistent) ways: both from a direct 
measurement of periods between threshold crossing, and from the Fourier power spectrum of the temporal signal 



The direct measurements of the T^pp from the time trace of p 2 {t) have been performed by defining two thresholds 



A. 



Distributions of the field intensities 



B. Average values of the LFF times 



p 2 (t). 
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FIG. 12: Field intensities distributions P(p 2 ) for the numerical data obtained by the integration of a noisy LK eqs. filtered at 
200 MHz. The data refer from left to rigth to p = 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, and 0.99 for (a) a = 3.3 
and (b) a = 4.0. Both the sets of data correspond to k — 0.35 and noise variance R — 3 x 10 -3 . 
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FIG. 13: Field intensities distributions P{p 2 ) for the numerical data obtained by the integration of a deterministic LK eqs. 
filtered at 200 MHz. The data refer from left to rigth to p = 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, and 0.99 for a = 5.0 
and k — 0.35. 



Ti < T2 and by idcntifing two consecutive time crossing of Ti, provided that in the intermediate time the signal 
has overcome the threshold Y2 at least once. The thresholds have been defined as Y\ =< p 2 > —2 * STD and 
T2 =< p 2 > +STD/2, where < ■ > and STD indicate the average and the standard deviation of the signal itself. 

The measurement of < T^pp > in terms of the power spectrum has been obtained by considering the power spectrum 
S(u>) of p 2 (t) and by evaluating the position u>m of the peak with the highest frequency, then < T^pp >= 2it/ujm- 
As already mentioned the two estimations are generally in very good agreement. 

In Fig. E| the average times < Tppp > are reported for two different sets of experimental measurements as a 
function of the pump parameter p, and compared with numerical data. In Fig. 1141 (a) are reported the experimental 
findings already shown in |lfil |. the estimation of Tlff have been performed both by direct inspection of the signal 
and via the first zero of the autocorrelation function (this second method corresponds to an evaluation from the 
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FIG. 14: (Color online) Average LFF times < Tlff > as a function of the pump parameter /i — fi red : black filled circles refer 
to experimental data, while the other symbols to the results obtained from the integration of the LK eqs. In the two figures 
are reported two different sets of experimental measures and the associated numerical data refer to k = 0.25 (a) and k = 0.35 
(b). In particular, red empty triangles correspond to the evolution of the noisy LK at a — 3.3 (with (a) /i red = 0.914 and 
(b) fi red =0.880) and blue empty squares to a = 4.0 (with fi red = 0.913 and 0.878, resp. for (a) and (b)). In both cases 
R = 3 x 10" 3 . The green stars denote the data of the deterministic LK eqs for a — 5.0 (in this case /i' e = 0.915 and 0.882, 
resp. for (a) and (b)). For the experimental measures fj, red — 0.916 in (a) and 0.875 in (b). The vertical dashed lines indicates 
the position of the solitary threshold for the experimental data, while the dash-dotted lines represent the decay c/(/x — fJ. red ), 
with c = 6.2 ns and 11 ns in (a) and (b), respectively. fi red is defined as the ratio /Ith- 



Fourier power spectrum). The numerical data have been obtained with the 2 methods outlined above for k = 0.25 for 
both noisy and deterministic LK eqs. In Fig. 1141 (b) a new set of experimental data is reported and compared with 
simulation results for k — 0.35, in this case all the data have been obtained by the method of the thresholds. From 
the figures it is clear that a reasonably good agreement between experimental and numerical data is observed for the 
deterministic case only for a = 5 (results for smaller a-values, obtained during the transient preceeding the stable 
phase, are not shown but they exhibits a worst agreement with experimental findings) and for the noisy dynamics for 
a = 4.0. 
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FIG. 15: (Color online) Standard deviation of the LFF times V as a function of the pump parameter /i — /i red : the symbols 
are the same as those reported in Fig. 1141 fbL The dash-dotted line indicates the power-law decay l/(/i — //" ed ) 3//2 . 
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A more detailed analysis can be obtained by considering not only the average values of the LFF times, but also the 
associated standard deviation V. This quantity reported in Fig. 1151 exibits a clear decrease with /i by approaching 
the solitary treshold, indicating a modification of the observed dynamics that tends to be more "regular" . Also in 
this case the comparison of experimental and numerical data suggests that the best agreement is again attained with 
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the noisy dynamics at a — 4.0. 

At this stage of the comparison we can sketch some preliminary conclusions: the LK eqs are able to reproduce 
reasonably well the experimental data for the VCSEL below the solitary threshold both in the deterministic case and 
in the noisy situation. However in the deterministic case a quite large value of the linewidth enhancement factor (with 
respect to the experimentally measured one) is required. A more detailed comparison will be possible by considering 
the PDFs of the T LFF . 



C. Distributions of the LFF times 



In this sub-section we will examine the whole distribution of the Tlff in more details. Considering the experimental 
data, we observe that all the measured PDFs obtained for different pump currents reveal an exponential-like tail at 
lo ng t imes and a rapid drop at short times (as shown in Fig. I16|) . These results are in agreement with those reported 
in 1 18 ] for a single-transverse-mode semiconductor laser in proximity of 1th ■ 
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FIG. 16: Probability density distributions of the Tlff- Solid lines refer to experimental data, the dashed ones to the Inverse 
Gaussian Distribution 1181 with the average and standard deviation corresponding to the experimental ones, (a) / = 2.56 mA, 
(b) I = 2.64 mA, and (c) J = 2.70 mA 



The typical dynamics corresponding to a LFF can be summarized as follows: a sudden drop of the intensity is 
followed by a steady increase of p 2 , associated to fluctuations of the intensity, until a certain threshold is reached 
and the intensity is reset to its initial value and restart with the same "sysiphus cycle" 0. This behaviour and the 
observed shapes of the PDFs suggest that the dynamics of the intensities can be modelized in terms of a Brownian 
motion plus drift. In other words, by denoting with x(t) the intensity, an effective equation of the following type can 
be written to reproduce its dynamical behaviour: 



(16) 



with initial condition x(0) = xq, where £(i) is a Gaussian noise term with zero average and unitary variance, rj repre- 
sents the drift, a is the noise strength. Within this framework the average first passage time to reach a fixed threshold 
T is simply given by r = (T — x )/r] , while the corresponding standard deviation is V = [\J (T — xo)cr]/r] 3 / 2 [I^. A res- 
onable assumption would be that ij is directly proportional to the pump parameter (/i— p red ), (where /i red = I^ d /Ith is 
the rescaled pump current value at the reduced threshold) and by further assuming that the threshold T is independent 
on the pump current this would imply that 



< T LFF >= 



- M red ) 



and Vlff = 



red\3/2 



(17) 



these dependences are indeed quite well verified for the experimental data above the solitary threshold as shown in 
Figs HI and Hoi 

For the simple model introduced by eq. (|16fl . the PDF of the first passage times is the so-called Inverse Gaussian 
Distribution [43| : 



P(T) 



^/27T 7 T 3 



-(T-t) 2 /(2jT) 



(18) 



where 7 = V 2 jr. A comparison of this expression with the experimentally measured P{Tlff) is reported in Fig. 
Hoi The good agreement suggests that the "sysiphus cycles" can be due to few elementary ingredients: a stochastic 
motion subjected to a drift plus a reset mechanism once the intensity has overcome a certain threshold. 
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A way of rewriting the distribution 118(1 in a more compact form as a function of only one parameter, the so-called 
coefficient of variation 5 — V/t (i.e., the ratio between the standard deviation V and the mean r), is to rescale the 
time as z — (T — t)/V and the PDFs as g(z) = V * P(T). This procedure leads to the following expression: 

g(*)= /o = = e-* 2 W+V . (19) 
■\j2ir{5z + l) 3 

It is clear that all the PDFs will coincide, once rescaled in this way, if the coefficient of variation 6 has the same value 
for all the considered pump currents. However this is not the case and indeed we measured values of S in the range 
[0.28 : 0.66] for I < I t h, nonetheless if we report in a single graph all these curves the overall matching is very good, 
as shown in Fig. 1171 




FIG. 17: Rescaled probability density distributions g(z) = V * P{Tlff) as a function of z = {Tlff— < Tlff >)/V. The 
curves refer to I = 2.54 mA, 2.56 mA, 2.64 mA, 2.70 mA, 2.75 mA, 2.80 mA, and 2.90 mA. 

Let us finally compare these distributions g(z) with the corresponding ones obtained from direct simulations of the 
LK eqs.. As one can see from Fig. El the agreement is good for the data obtained from the simulation of the noisy 
LK eqs at a — 4.0, while it is worse for the deterministic LK eqs at a = 5.0. 



VIII. CONCLUSIONS 



We presented a detailed experimental and numerical study of a semiconductor laser with optical feedback. The 
choice of a Vertical Cavity Laser pumped close to its threshold, together with a polarized optical feedback, assures 
a great control over the possibility of lasing action of other orders longitudinal and/or transverse modes than the 
fundamental one, and of the activation of the other polarization. In such a way, the description of the system using 
the Lang-Kobayashi model is well justified and it allows for a meaningful comparison with the experimental data. 
The analysis has been performed with particular regard to the LFF regime, where the model has been numerically 
integrated using parameters carefully measured in the laser sample used for the measurements. 

The comparison of the the measurements carried out in the VCSEL with polarized optical feedback with the 
predictions of the deterministic LK model suggest that in the examined range of parameters the dynamics of the 
model is characterized by a chaotic transient leading to stable ECMs with high gain. The transient duration increases 
(and possibly diverges) with increasing values of the rescaled pump current /i and of the linewidth enhancement factor 
a. We have not found evidence of periodic or quasi-periodic asymptotic attractors, as instead reported in [||, this can 
be due to the a-range examined in the present paper (namely 2.4 < a < 5.5), since these solutions become relevant 
for the dynamics only for a > 5 (as stated in [6j). 

However, a stationary LFF dynamics with characteristics similar to those measured experimentally, can be obtained 
for realistic values of the a parameter (namely, a ~ 3 — 4) only via the introduction of an additive noise term in the 
LK equations. The role of noise in determining the statistics and the nature of the dropout events has been previously 
examined in 0, Q , but in the present paper we have clarified that the LFF dynamics can be interpreted at a first 
level of approximation as a biased Brownian motion towards a threshold with a reset mechanism. 
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FIG. 18: (Color line) Rescaled probability density distributions g(z) = V *P(Tlff) as a function of z = (Tlff— < Tlff >)/V. 
The thick solid curve refer to experimental data for / = 2.75 mA (corresponding to /i ~ 1), the other ones to numerical findings: 
(a) data obtained for the noisy LK eqs for a = 4.0, k = 0.35 with variance R = 3 x 10 -3 and corresponding to fi — 0.90 (red), 
0.94 (blue), 0.98 (green), and 1.02 (violet); (b) data for the the deterministic LK eqs. for a — 5.0, k — 0.35 corresponding to 
/i = 0.92 (red) , 0.94 (blue), 0.96 (green), 0.98 (violet), and 1.00 (orange). 
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